Identifying the impacts of extreme weather

A case study of Houston, Texas

Link to GitHub repository: https://github.com/aakriti-poudel-chhetri/impacts-of-extreme-weather

Background

From February 13 to 17, 2021, Winter Storm Uri brought unusually cold temperatures, heavy snowfall, ice and freezing rain across Texas. Many places experienced their coldest conditions in more than three decades. The state’s power infrastructure was not designed for such extreme cold, therefore numerous power plants and grid components froze or malfunctioned. At the same time, electricity demand surged as residents attempted to heat their homes, placing severe pressure on the power system. To avoid a complete grid failure, authorities implemented rolling blackouts that were intended to last about an hour, but in reality many households lost electricity for several hours or even longer. An estimated 4.5 million homes experienced outages. These disruptions triggered additional problems, including widespread water and food shortages, burst pipes, and a lack of indoor heating. At least 246 deaths were linked directly or indirectly to the storm. Economic losses exceeded $195 billion, making Winter Storm Uri the expensive natural disaster in Texas history.

Understanding how communities were affected is crucial for preparing for future extreme weather events. It can help identify which neighborhoods are most vulnerable and guide decisions about where to focus recovery and emergency resources.

Thus, this study examines how Winter Storm Uri affected communities across the Houston metropolitan area by estimating the number of homes that lost power and assessing whether lower-income neighborhoods faced a disproportionate share of the impacts. Using satellite-based nighttime lights data from the VIIRS VNP46A1 product, we identify areas where electricity was likely lost by comparing light intensity before and after the storm. These blackout areas are then combined with OpenStreetMap building data to estimate the number of affected homes and linked with U.S. Census Bureau data to evaluate socioeconomic differences—such as variations in median household income—between impacted and non-impacted census tracts.

Research Question

Understanding how power outages unfolded across Houston during Winter Storm Uri is important for evaluating both the extent of the blackout and the communities most heavily affected. Using satellite-based detection of outage areas alongside socioeconomic data, this study examines whether the storm’s impacts were distributed unevenly across neighborhoods. The main research question guiding this analysis is:

“How were power outages during Winter Storm Uri distributed across the Houston metropolitan area, and did these outages disproportionately affect lower-income neighborhoods?”

This study also addresses the following sub-questions:

• How many homes likely lost electricity during the storm?
• Which census tracts were most affected by the blackout?
• How does blackout exposure vary with median household income across census tracts?

Data

Data source

For this analysis, data were extracted from the multiple sources.

Night lights

• The NASA’s Worldview is explored for extracting the data around the day of the storm. There are several days with too much cloud cover to be useful, but 2021-02-07 and 2021-02-16 provides two clear, contrasting images to visualize the extent of the power outage in Texas.

VIIRS data is distributed through NASA’s Level-1 and Atmospheric Archive & Distribution System Distributed Active Archive Center (LAADS DAAC). Many NASA Earth data products are distributed in 10x10 degree tiles in sinusoidal equal-area projection. Tiles are identified by their horizontal and vertical position in the grid. Houston lies on the border of tiles h08v05 and h08v06. These two tiles per date were pre-downloaded and provided by the team.

Roads

• We used Geofabrik’s download sites to retrieve a shapefile of all highways in Texas and prepared a Geopackage (.gpkg file) containing just the subset of roads that intersect the Houston metropolitan area.

Houses

• The data were downloaded from Geofabrik and processed to create a GeoPackage containing only residential buildings within the Houston metropolitan area.

Socioeconomic

• The socioecoomic data were obtained from the U.S. Census Bureau’s American Community Survey for census tracts in 2019.

Data analysis process

To assess the impacts of Winter Storm Uri on Houston communities, a multi-step approach combining spatial and socioeconomic data was adopted. This allowed us to identify areas affected by power outages, estimate the number of homes impacted, and examine patterns across different socioeconomic groups. The major steps of the analysis are summarized below:

• Locate outrage areas: Using VIIRS VNP46A1 satellite-based nighttime lights data to identify areas likely experiencing power outages by comparing light intensity before and after the storm.
• Estimate affected homes: Overlay blackout zones with building footprint data from OpenStreetMap to estimate the number of homes without electricity.
• Socioeconomic analysis: Link spatial data with U.S. Census Bureau information to examine differences in median household income between impacted and non-impacted census tracts.
• Combine assessment: Combine spatial and socioeconomic insights to provide a comprehensive understanding of how the storm affected local communities.

Set up

Let’s start by importing all the necessary libraries.

library(tidyverse)
library(sf)
library(tmap)
library(terra)
library(stars)
library(kableExtra)

First, we read our datasets.

# Read night lights data
tile05_07 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021038.h08v05.001.2021039064328.tif'))
tile06_07 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021038.h08v06.001.2021039064329.tif'))
tile05_16 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather','data', 'VNP46A1', 'VNP46A1.A2021047.h08v05.001.2021048091106.tif'))
tile06_16 <- read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021047.h08v06.001.2021048091105.tif'))


# Read roads data
highways <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'gis_osm_roads_free_1.gpkg'),
                 query = "SELECT * FROM gis_osm_roads_free_1 WHERE fclass = 'motorway'") %>%
  st_transform(crs = 'EPSG:3083')


# Read houses data
houses <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data',
                          'gis_osm_buildings_a_free_1.gpkg'),
  query = "SELECT *
    FROM gis_osm_buildings_a_free_1
    WHERE (type IS NULL AND name IS NULL)
      OR type IN ('residential', 'apartments', 'house', 'static_caravan', 'detached')") %>% 
    st_transform('EPSG:3083')


# Explore the contents of the geodatabase
socioeconomic <- st_layers(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'))

# Extract the census tract information
census_tract <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'), layer = 'ACS_2019_5YR_TRACT_48_TEXAS') %>% 
  st_transform('EPSG:3083')

# Extract the income information
income <- st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'), layer = 'X19_INCOME')

Night light intensity before & after the winter storm

We will merge the raster tiles and crop it to the Houston bounding box to focus the analysis on the study area. It will ensure only valid night-light data within Houston were used for assessing power outages.

# Merge two raster data of February 7
nightlight_07 <- st_mosaic(tile05_07, tile06_07)

# Merge two raster data of February 16
nightlight_16 <- st_mosaic(tile05_16, tile06_16)

# Create Houston bounding box
houston_bbox <- st_bbox(c(xmin = -96.5, xmax = -94.5, ymin = 29, ymax = 30.5),
                        crs = 'EPSG:4326')

# Crop all raster data to Houston area
nightlight_07_crop <- st_crop(nightlight_07, houston_bbox)
nightlight_07_crop[(nightlight_07_crop > 1000) | (nightlight_07_crop <= 0)] <- NA

nightlight_16_crop <- st_crop(nightlight_16, houston_bbox)
nightlight_16_crop[(nightlight_16_crop > 1000) | (nightlight_16_crop <= 0)] <- NA

Here, we visualize the comparison between night-light intensities in Houston before and after the February 2021 winter storm. By mapping the February 7 (pre-storm) and February 16 (post-storm) night-light data side by side, the analysis shows areas where lighting and likely power availability changed, helping to identify regions affected by the blackout.

tmap_options(component.autoscale = FALSE)

# Map for February 7 (before storm)
beforestorm <- tm_shape(nightlight_07_crop) +
  tm_raster(col.scale = tm_scale_continuous(values = 'inferno'),
            col.legend = tm_legend(show = TRUE, 
                                   title = 'Night light intensity before storm',
                                   title.size = 1,
                                   orientation = 'landscape',
                                   width = 35,
                                   height = 3,
                                   legend.just = c(0, 1)))+
  tm_title('Houston Night Lights\nFebruary 07, 2021 (Before storm)',
           size = 1.2, fontface = 'bold') +
  tm_compass(position = c('left', 'bottom'),
             type = 'arrow',
             size = 1.5,
             text.size = 0.7,
             color.dark = 'grey50',
             text.color = 'white') +  
  tm_scalebar(position = c('left', 'bottom'),
              text.size = 0.8,
              color.dark = 'grey50',
              text.color = 'white') +  
  tm_layout(bg.color = 'white',           
            outer.bg.color = 'white',      
            frame = TRUE,
            legend.show = TRUE) +
  tm_graticules(x.show = TRUE, y.show = FALSE, lwd = 0.2, col = "white", alpha = 0.3, labels.size = 0.5)
  
  
  #tm_graticules(lwd = 0.2, col = "white", alpha = 0.3, labels.size = 0.5)


# Map for February 16 (after storm)
afterstorm <- tm_shape(nightlight_16_crop) +
  tm_raster(col.scale = tm_scale_continuous(values = 'inferno'),
            col.legend = tm_legend(show = TRUE, 
                                   title = 'Night light intensity after storm',
                                   title.size = 1,
                                   orientation = 'landscape',
                                   width = 35,
                                   height = 3)) +
  tm_title('Houston Night Lights\nFebruary 16, 2021 (After storm)',
           size = 1.2, fontface = 'bold') +
  tm_compass(position = c('left', 'bottom'),
             type = 'arrow',
             size = 2.5,
             text.size = 0.7,
             color.dark = 'grey50',
             text.color = 'white') +  
  tm_scalebar(position = c('left', 'bottom'),
              text.size = 0.8,
              color.dark = 'grey50',
              text.color = 'white') +  
  tm_layout(bg.color = 'white',           
            outer.bg.color = 'white',      
            frame = TRUE,
            legend.show = TRUE) +
  tm_graticules(x.show = TRUE, y.show = FALSE, lwd = 0.2, col = "white", alpha = 0.3, labels.size = 0.5)

tmap_arrange(beforestorm, afterstorm, ncol = 2)

Figure 1: Comparison of Houston night light intensity before and after the winter storm, February 2021

Identify the area with power outage

We will now create a blackout mask by identifying areas of Houston where night light intensity decreased significantly after the storm, indicating potential power outages. By calculating the difference between before and after night light, cropping to the study area, filtering out minor changes, and vectorizing the result, the mask provides a spatial representation of blackout affected regions for further analysis.

# Calculate the difference in light intensity
nightlight_diff <- nightlight_07 - nightlight_16

# Crop and reclassify the raster data 
mask <- st_crop(nightlight_diff, houston_bbox)
mask[mask < 200] <- NA

# Vectorize the blackout mask
blackout <- st_as_sf(mask,
                     as_points = FALSE,
                     merge = TRUE) %>%
  st_make_valid() %>%
  st_transform(crs = 'EPSG:3083')

Exclude highways from the power outage area

In this step, highway geometries are combined and buffered to exclude areas where changes in night light intensity could be due to traffic rather than power outages. By removing blackout areas within 200 meters of highways, the analysis isolates regions where lighting loss more likely reflects true residential or commercial power outages.

# Combine all highway geometries into one
highways_union <- st_union(highways)

# Create 200m buffer around all highways
highways_buffer <- st_buffer(highways_union, dist = 200)

# Find areas that experienced blackouts further than 200m from highways
blackout_far <- st_difference(blackout, highways_buffer)

We visualize the spatial extent of power outages across Houston by mapping the blackout areas. The map highlights the neighborhoods affected, providing a clear view of the geographic distribution of the electricity disruption during the storm.

# Visualize the blackout in Houston
tm_basemap('OpenStreetMap') +
  tm_shape(blackout_far) +
    tm_polygons(fill = 'ivory', fill_alpha = 0.8, col = 'black') +
    tm_title('Blackout areas in Houston', fontface = 'bold', size = 1) +
    tm_layout(legend.show = TRUE) +
    tm_scalebar(position = c('left', 'bottom')) +
    tm_compass(position = c('right', 'top')) +
  tm_graticules(lwd = 0.2, alpha = 0.3, labels.size = 0.5) +
  tm_add_legend(labels = c("Blackout Areas"),
              fill = c('black'),      
              type = 'polygons',
              position = c('left', 'top'))

Figure 2: The figure shows the spatial distribution of power outages in Houston, highlighting areas impacted by blackout.

Identify number of homes impacted by power outage

We estimate the number of homes in Houston that lost power during the February 2021 winter storms. By identifying which building footprints intersect with the blackout mask, the analysis isolates the homes likely affected by outages. Converting buildings to centroids allows for precise spatial calculations and comparisons between impacted and non-impacted houses, resulting in an estimated total of approximately 157,970 homes without power.

# Keep buildings that intersect with blackout areas
houses_blackout <- st_intersects(houses, blackout_far)

# Check if each building is in a blackout area
in_blackout <- lengths(houses_blackout) > 0

# Filter buildings in blackout areas
blackout_houses <- houses[in_blackout, ]

# Number of impacted buildings
n_blackout_houses <- nrow(blackout_houses)

# Convert all buildings to centroids
houses_points <- houses %>% 
  st_centroid()

# Houses impacted by blackouts
houses_impacted <- blackout_houses %>% 
  st_centroid()

# Houses not impacted by blackouts
houses_not_impacted <- houses_points %>% 
  filter(!osm_id %in% blackout_houses$osm_id)

Let’s create a table to show the number of houses affected by the power outage, making it easier to visualize the impact.

# Table of homes that experienced blackouts
houses_df <- houses_impacted |> 
  group_by(type) |> 
  summarise(count = n(),
            percent = round((n() / nrow(houses_impacted)) * 100, digits = 2)) |>
  st_drop_geometry()

# Add total row
houses_df <- houses_df |> 
  add_row(type = 'Total', count = sum(houses_df$count))

# Change type NA to character string
houses_df$type[6] <- "NA"

# Display kable table
options(knitr.kable.NA = "")
kbl(houses_df,
    table.attr = 'class="custom-table"',
    col.names = c(" Types", "Count", "Percentage (%)"),
    align = 'c')

Types	Count	Percentage (%)
apartments	1136	0.72
detached	353	0.22
house	19760	12.51
residential	1395	0.88
static_caravan	80	0.05
NA	135243	85.61
Total	157967

Plot map for the houses that were impacted by the power outage.

# Create the map
tm_basemap('OpenStreetMap') +
  tm_shape(houses_not_impacted) +
    tm_dots(fill = 'midnightblue', size = 0.02, fill_alpha = 0.8) +
  tm_shape(houses_impacted) +
    tm_dots(fill = 'firebrick', size = 0.02, fill_alpha = 0.8) + 
  tm_title('Houses in Houston that lost power',
           position = tm_pos_out('center', 'top'), 
           fontface = 'bold', size = 1.2) +
  tm_layout(legend.outside = TRUE,
            legend.outside.position = 'right',
            legend.title.size = 1,
            legend.text.size = 0.8) +
  tm_add_legend(labels = c('Impacted houses', 'Not impacted houses'),
                fill = c('firebrick', 'midnightblue'),
                col = c('firebrick', 'midnightblue'),
                alpha = c(0.2, 0.7, 0.5),
                title = 'Legend') +
  tm_scalebar(position = c('left', 'bottom'),
              text.size = 0.6) +
  tm_compass(position = c('right', 'bottom'),
             type = '4star',
             size = 2)

Figure 3: This map visualizes the spatial distribution of homes in Houston that lost power during the February 2021 winter storms. Impacted houses are highlighted in red, while unaffected houses are shown in blue, allowing for a clear comparison of affected versus unaffected areas across the city. The map provides a detailed spatial perspective on the extent of the blackout and helps identify which neighborhoods were most impacted.

Identify census tracts impacted by power outage

The following process links census tract boundaries with median household income data to analyze the socioeconomic characteristics of neighborhoods affected by the blackout. By joining the income dataset to the census tract geometries, the analysis enables comparisons between impacted and non-impacted areas and supports investigation of whether power outages disproportionately affected higher- or lower-income communities.

# Rename the column in income dataset to match with the census tract
income_renamed <- income %>% 
  rename("GEOID_Data" = "GEOID")

# Left join the two datasets by the geometry
census_income <- census_tract %>% left_join(income_renamed, by = 'GEOID_Data')

Let’s verify CRS. Verifying CRS alignment is essential for accurate spatial analysis, as mismatched projections can lead to incorrect geometry intersections, area calculations, or mapping errors.

# Check CRS of the dataset
if(st_crs(census_income) == st_crs(houston_bbox)){
print("Coordinate reference systems match!")
} else {
warning("Update coordinate reference systems to match!")
}

Warning: Update coordinate reference systems to match!

By aligning the coordinate systems and cropping to the Houston study area, census tracts and blackout affected houses are prepared for spatial comparison. Intersections between the two datasets identify which tracts experienced power outages, and a new column flags these impacted areas. Filtering to only the affected tracts enables a targeted analysis and visualization of how the blackout influenced neighborhoods across the city.

# Transform the CRS and crop to houston bounding box
census_income <- st_transform(census_income, crs = 'epsg:4326') 

# Crop census income to Houston
census_income_crop <- census_income %>% 
  st_crop(houston_bbox)

# Transform and clean blackout houses
blackout_houses_transformed <- blackout_houses %>% 
  st_transform(crs = 'epsg:4326')

# Identify the blackout areas
census_blackout <- st_intersects(census_income_crop, blackout_houses_transformed)

# Create the column
census_blackout_col <- census_income_crop %>% 
  mutate(blackout = lengths(census_blackout) > 0)

# Select only the values that are true for plotting
census_blackout_col_filtered <- census_blackout_col[census_blackout_col$blackout, ]

Visualize the census tracts in Houston that experienced power outages on a map.

# Map the census tracts that lost power
tm_basemap('CartoDB.VoyagerNoLabels') +
  tm_shape(census_blackout_col_filtered) +
  tm_polygons(fill = 'blackout',
              fill.scale = tm_scale_categorical(
                values = c('TRUE' = 'firebrick'),
                labels = c('TRUE' = 'Blackout')),
              fill.legend = tm_legend(title = 'Blackout status'),
              col = 'grey50',
              lwd = 0.2) +
  tm_title('Census tracts in Houston that lost power',
           position = tm_pos_out('center', 'top'),
           fontface = 'bold',
           size = 1) +
  tm_layout(legend.outside = TRUE,
            legend.outside.position = 'right',
            legend.bg.color = 'white',
            legend.bg.alpha = 0.8) +
  tm_compass(position = c('right', 'bottom'),
             type = 'arrow',
             size = 1) +
  tm_scalebar(position = c('left', 'top'),
              text.size = 0.8)

Figure 4: The following map highlights Houston census tracts that experienced power outages during the February 2021 winter storms. Affected tracts are shown in red, providing a clear visualization of the spatial distribution of blackouts across the city. This allows for easy identification of neighborhoods impacted by the storm and supports further analysis of socioeconomic patterns in relation to power loss.

We compare the distributions of median household income between Houston census tracts that experienced blackouts and those that did not during the February 2021 winter storms. By comparing the differences in income, it allows for assessment of whether power outages disproportionately affected higher or lower income neighborhoods, and provides insight into potential socioeconomic disparities in storm impacted areas.

ggplot(census_blackout_col %>% filter(!is.na(B19013e1)), # Remove NA values
                                      aes( x = blackout, y = B19013e1, fill = blackout)) +
  geom_boxplot(alpha = 0.8) +
  scale_x_discrete(labels = c('TRUE' = 'Blackout',
                              'FALSE' = 'No blackout')) +
  scale_fill_manual(values = c('TRUE' = 'firebrick',
                               'FALSE' = 'midnightblue'),
                    labels = c('Blackout', 'No blackout')) +
  labs(title = 'Median household income for census tracts',
       subtitle = 'Blackout status in Houston census tracts',
       x = 'Blackout status',
       y = 'Median household income in $',
       fill = 'Status') +
  scale_y_continuous(labels = scales::dollar_format()) +
  theme_minimal() +
  theme(plot.title = element_text(face = 'bold', size = 16),
        plot.subtitle = element_text(size = 14),
        legend.position = 'bottom')

Figure 5: This plot compares the distributions of median household income between census tracts that experienced blackouts and those that did not during the severe winter storms of February 2021 in Houston, Texas.

Conclusion

This analysis used satellite nighttime light data to identify which parts of Houston lost power during Winter Storm Uri and how these outages related to neighborhood income levels. Surprisingly, the results showed that higher income census tracts appeared slightly more affected, even though past research, including Lee et al. (2022), found that low income and minority communities experienced more severe and prolonged outages. This difference likely comes from limitations in our dataset. Wealthier neighborhoods often have more outdoor lighting, which can make blackout detection less accurate. Census tract data also hides smaller neighborhood differences, and highway lighting or cloud cover may have caused some areas to be misclassified. In addition, assuming all buildings were residential and using a 200-meter buffer around highways may have excluded areas that actually lost power. Our analysis estimated that about 157,967 residential buildings lost electricity. However, because only two days of satellite data were used, the results may not capture the full timeline of outages. A longer range of dates would make it easier to see where and when blackouts occurred, and provide a clearer picture of storm impacts.

Next steps

While blackout exposure appeared similar across income levels, the ability to cope with and recover from outages is not equal. Low income communities have fewer resources to stay warm, relocate or rebound after a long outage. This difference between exposure and impact is important for understanding who is most vulnerable during extreme weather events. Future work should use nighttime light data from more days to better track outage duration and recovery. Adding other socioeconomic information, such as age, housing quality, or access to transportation, can improve our understanding of community vulnerability. Combining multiple sources like utility outage records, mobile phone mobility data and satellite imagery would create a more complete picture of who was affected and how. Examining infrastructure characteristics, such as the age of power lines or the placement of substations, may also help explain differences in recovery across neighborhoods. Together, these steps can help identify communities at higher risk and support more effective and equitable preparation for future extreme weather events.

References